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ABSTRACT 



Context. For most astronomical objects, radiation is the only probe of their physical properties. Therefore, it is important to have 
the most elaborate theoretical tool to interpret observed spectra or images, thus providing invaluable information to build theoretical 
models of the physical nature, the structure, and the evolution of the studied objects. 

Aims. We present IRIS, a new generic three-dimensional (3D) spectral radiative transfer code that generates synthetic spectra, or 
images. It can be used as a diagnostic tool for comparison with astrophysical observations or laboratory astrophysics experiments. 
Methods. We have developed a 3D short-characteristic solver that works with a 3D nonuniform Cartesian grid. We have implemented 
a piecewise cubic, locally monotonic, interpolation technique that dramatically reduces the numerical diffusion effect. The code 
takes into account the velocity gradient effect resulting in gradual Doppler shifts of photon frequencies and subsequent alterations of 
spectral line profiles. It can also handle periodic boundary conditions. This first version of the code assumes Local Thermodynamic 
Equilibrium (LTE) and no scattering. The opacities and source functions are specified by the user. In the near future, the capabilities 
of IRIS will be extended to allow for non-LTE and scattering modeling. 

Results. IRIS has been validated through a number of tests. We provide the results for the most relevant ones, in particular a searchlight 
beam test, a comparison with a ID plane-parallel model, and a test of the velocity gradient effect. 

Conclusions. IRIS is a generic code to address a wide variety of astrophysical issues applied to different objects or structures, such 
as accretion shocks, jets in young stellar objects, stellar atmospheres, exoplanet atmospheres, accretion disks, rotating stellar winds, 
cosmological structures. It can also be applied to model laboratory astrophysics experiments, such as radiative shocks produced with 
high power lasers. 

Key words, methods: numerical - radiative transfer 



1. Introduction 

Essentially all objects in the Universe have a complicated, 
three-dimensional, dynamic structure. Understandably, through- 
out most of the history of astronomy the observed objects had 
to be modeled using significant restrictions and approximations. 
Among them, a set of assumptions about the global geometry 
of an object always played a pivotal role. For instance, stellar 
atmospheres were typically treated assuming a plane-parallel, 
horizontally-homogeneous geometry, which simplified the prob- 
lem to one spatial dimension. Similarly, stellar winds, novae, 
planetary nebulae, and other extended sources, were modeled us- 
ing an assumption of spherical symmetry, which again renders it 
a one-dimensional problem. 

In the last two decades, there was much activity de- 
voted to extend the traditional modeling techniques to treat 
structures that are truly 3-dimensional (3D). In fact, detailed 
3D hydrodynamic simulations have now become essentially 
routine. A non-exhaustive list of 3D astrophysical (radi- 
ation) (magneto)hydrodynarmc^_(RHD 2 _MHD_or RMHD) 
codes includes ZEUS dStone & Normanl Il992allbl). ENZO 
(iBrvan & Normanlll997l: INorman & Brvanlll999t lO'Shea et al] 
120041: INorman et alJ 120071) RAMSES (iTevssieri 120021) 
HERACLES dGonzalez et all l2007h . PLUTO dMignone et alJ 



120071). ATHENA dStone et alJ 120081) and its radiation m odule 
dDavis etal.ll2012l) . AstroBEAR dCunningham et alJl2009h . 



However, since most of the information about an astronom- 
ical object is conveyed to a distant observer by its radiation, de- 
tailed 3D hydrodynamic models have to be accompanied by ade- 
quate 3D radiation transfer solutions that provide a spectroscopic 
diagnostic information. Compared to 3D hydrodynamic models, 
the 3D radiation transfer solvers are much more complicated and 
difficult to treat numerically, because (i) there are many more 
quantities that describe a radiation field, as compared to the hy- 
drodynamics, due to the directional and spectral dependence of 
radiation; (ii) a long-range interaction between the radiation field 
and the plasma arises because of a typically much larger mean 
free path of a photon compared to a mean free path of massive 
particles. 

Although providing exact time-dependent, non-LTE radi- 
ation (magneto)hydrodynamic models of astronomical objects 
is generally viewed as a mighty goal, such a goal has not 
been fully achieved yet. However, a large progress was accom- 
plished in recent years, specifically by 3D RMHD codes for stel- 
lar atmospheres, such a s: the Copenhagen-Oslo Stagger Codes 
dNordlund & Galsgaardl Il995b iG alsgaard & Nordlund 1996; 
| Hansteenl2004l:lHansteen et al.l2007l) CQ 5 BOL D dFrevtag et al 
l2002h . MURaM dVogler et all 120051) . Bifrost dGudiksen et al 
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1 d_ 

c dt 



+ n • V I(r, n, v, t) 



1201 ll) . All of these codes solve the 3D radiative transfer by iMihalas & Mihalasl [l984) 
approximating the monochromatic opacities with a s mall num- 
ber o f mean opacities (the opacity binning scheme of Nordlund 
Il982l) . With typically four bins, they can model the frequency- 
integrated radiative energy losses and gains quite well. 

There are also several recently developed computer 
codes specifically designed to provide 3D radiative trans- 
fer solvers, typically for a post-processing of 3D hy- 
drodynamic simulations. The most widely used among 



rj(r, n, v, t) - x( r , n , v > K r < n,v,t) , (1) 



1994 



1995 



them are the f ol lowing: MUGA (lAuer & Paletou 
Auer etal 1 Il994t If ruiillo Bueno & Fabiani Bendicho 
Fabiani Bendi cho et al l 1 19971). MULTI3D dBotnenl 
Leenaarts & Carlssonl 120091 : iLeenaarts et al. 120091). which is 



1997 



a generalizati on of the ID co d e MU LTI dCarlssonl Il986h . 
the RH code dUitenbroekl 1200 ll 120061). a 3D version of the 
atmospheric code PHOENIX dHauschildt & Baronl l2006t 
iBaron & Hauschildtl 120071 iHauschildt & Baronl 120081). ASSeT 
dKoesterke etal.1 120081: iKoesterkd 120091) . SCATE (iHaveketall 
|201j. They differ in their intended range of applications, and 
in many details of the numerical techniques. They all use some 
variant of the short-characteristics scheme (see Section [3), but 
differ in details of what is assumed about the behavior of state 
parameters between the grid points, and how the necessary 
interpolations are being performed. For details, the reader is 
referred to the above mentioned pa pers, and several excellent 
review papers dCarlssonll2008ll2009l) . 

In this paper, we describe our variant of the three- 
dimensional radiative transfer solver, named IRIS. Analogously 
to most of the above mentioned techniques, it uses the short- 
characteristics scheme. Our solver differs from the previous ones 
in several respects: it is formulated only in the Cartesian co- 
ordinate system, but within this system it is formulated in the 
most universal way that allows for several variants of the spa- 
tial integration and interpolation to be tested and used depend- 
ing on the actual application; it carefully treats subgriding to 
allow for line radiation transport in the presence of arbitrary 
(non-relativistic) velocity fields. Although the current version 
assumes Local Thermodynamic Equilibrium (LTE) and no scat- 
tering, the code is prepared to be relatively straightforwardly ex- 
tended to treat scattering and departures from LTE, the so-called 
non-LTE (or NLTE) situations. This will be reported in future 
papers. 

The paper is organized as follows. In Section [2] we outline 
the main physical assumptions of the model solved by IRIS. In 
Section [3] we summarize the short-characteristics method, and 
its application in a 3D Cartesian grid. In Section [4] we describe 
the piecewise, cubic, locally monotonic, interpolation technique 
that we use. Then, in Section [5] we describe our implementa- 
tion in IRIS of the techniques presented in Sections [3] and |4] 
We also explain the procedure to handle general macroscopic 
velocity fields. The calculation of the moments of the radiation 
field and the related angular integration methods are the subjects 
of Section [6] We detail in Section [7] the method that we em- 
ploy to treat media with periodic boundary conditions. Section[8] 
presents the results of relevant tests. We conclude our paper in 
Section [9] with a summary of the features of IRIS, and with an 
outline of possible extensions to the code. 



2. Radiative Transfer: Basic Definitions and 
Assumptions 

The general unpolarized r adiative transfer equation (RTE) in 
the observer's frame reads dChandrasekharl 1 9501: iMihalasll 1 9781: 



where I(r, n, v, t) is the specific intensity of radiation at position 
r, propagating in the direction specified by the unit vector n, 
at frequency v, and time f; ^(r, n, v, t) is the absorption coeffi- 
cient, rj(r, n, v, t) is the emission coefficient, and c is the speed of 
light. Although the above quantities depend on time, the current 
version of IRIS solves the time-independent RTE, for a given 
hydrodynamics structure at a given instant, i.e., for a given snap- 
shot. In other words, we consider regimes in which the photon 
free-flight time is small compared to the fluid flow dynamical 
timescales, so that the radiation field gets fully stabilized before 
any change occurs in the flow dynamical properties. We assume 
that the typical flow velocities of the hydrodynamic structures 
are non reiativistic. 

As is customary, we introduce the source function, 

T](r, n, v, f) 



S (r, n, v, t) 



(2) 



To simplify the notations, we drop the explicit dependence of ra- 
diative quantities on the position, direction and time, and denote 
the dependence on frequency with a subscript v. The transfer 
equation ([TJ becomes 

j^=Xv{S v -I v ), (3) 

where s is the path length along the ray in the direction of prop- 
agation of the radiation. 

The choice of a numerical method depends on the purpose of 
the simulation. If one is interested only in computing emergent 
specific intensity from a 3D computational box, and if the opac- 
ity and the source function are fully known within the box (the 
so-called "formal solution"), one can select a number of rays 
(photon paths) that emerge from the box and cross the whole 
extent of the box, and then solve the radiative transfer equation 
along such rays. Obviously, a ray does not generally go through 
the original grid points. One can choose any discretization of 
the ray, and then obtain needed values of the opacity and source 
function at the discretized positions by a three-dimensional inter- 
polation. A more efficient way is to discretize the ray by taking 
its intersections with the individual planes defined by the grid 
points; in this case one deals with two-dimensional interpola- 
tions. 

The problem is thus reduced to a set of ID problems. Any 
method capable of solving the transfer equation along a ray 
can be used here; most popular in the astrophysica l applica- 
tions being the Feautrier metho d, its 2nd ord er variant ( Fe autrierl 



1964) or the 4-th order variant (|Aue 
Finite Element method (Castor et 



jen 

in 



1976); the Discontinuous 
1992); or the "ID short- 



characteristics" method dOlson & Kunaszll987l) . All these meth 
ods are called the "long-characteristics" scheme in some specific 
situations. However, the terminology is not used consistently; 
some studies use this term in a more restricted meaning, namely 
for using a ID short-characteristics scheme for solving the trans- 
port for a set of rays passing through the whole computational 
region. 

The situation is different when the source function is not 
known a priori within the computational box. The simplest situ- 
ation of this sort arises when the opacity is specified in the com- 
putational domain and one can even assume LTE, but the contin- 
uum scattering is taken into account. Another case is a general 
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non-LTE situation where line scattering also is important. The 
most complex situation is when dealing with the true radiation 
hydrodynamics, where one has to pass the values of radiation 
moments at all grid points to the hydrodynamics part of the code. 

Even in the simplest case of LTE with electron (Thomson) 
scattering, the total source function depends on the mean inten- 
sity of radiation, and therefore has to be determined iteratively. 
The essential point is, however, that even if we are interested 
only in the emergent radiation, we have to determine the total 
source function, and therefore the specific intensities, in all grid 
points to be able to interpolate it into the discretized positions 
along the rays. 

Although the long-characteristics scheme is useful for some 
purposes, such as allowi ng for an efficient paralle l scheme within 
domain decomposition (Heinemann et al. 2006), we chose the 
short-characteristics scheme. The classical setup of the long- 
characteristics method is to pass a ray through each grid point 
through the whole computational box, in which case the number 
of necessary operations is obtained as follows. Let us assume 
that the computational box is discretized using N points in each 
direction, so that there are A^ 3 grid points (and also « A^ 3 el- 
ementary cells). A long-characteristics scheme would consider 
Af 3 rays. Solving the transfer equation along one ray, for one di- 
rection and one frequency, would require O(N) operations, an 
exact value depending on how exactly are the necessary interpo- 
lations performed, assuming that they require 0(1) operations. 
For one direction and one frequency one thus needs 0(N 4 ) oper- 
ations. In contrast, the short-characteristics scheme solves the 
transfer equation within one cell only, and therefore the total 
number of operations is proportional to the number of cells, that 
is 0(N 3 ) operations. 

To avoid confusion, we mention that when using the long- 
characteristics scheme for obtaining an emergent intensity only, 
it requires also (9(A^ 3 ) operations, because we have A^ 2 rays; 
each going through one grid point on the boundary plane that 
faces an "observer", and a transfer solution along each ray re- 
quires O(N) operations. In this particular case, the scaling of 
the computer time is the same as for the short-characteristics 
scheme. In some cases, the use of long-characteristics scheme 
may even be more computationally advantageous than the short- 
characteristics, and therefore some codes (e.g., ASSeT, SCATE) 
use the short-characteristics scheme when dealing iteratively 
with scattering, and use the long-characteristics scheme to eval- 
uate the emergent intensity for the final converged model. 

We note that there are several possible numerical methods 
that propagate an information form one edge of the cell to an- 
other, whose numbers of operations also scale as the number of 
cells, 0(N 3 ). In astrophysics, the vast majority of approaches is 
based on the short-characteristics scheme, and this is what we 
adopt in this paper as well. In this paper we assume LTE and no 
scattering. In this case, the absorption coefficient Xv is a function 
of the local mass density p and the local temperature T, and the 
source function S v is equal to the local Planck function B V (T). 
With these two assumptions (LTE and no scattering), we could 
have adopted the long-characteristics method as well, but in view 
of our intended future development of the solver to treat scatter- 
ing and non-LTE situations, the short-characteristic scheme is 
clearly the appropriate choice. 



3. 3D Short-Characteristics 

In the conte xt of astrophys i cal rad iative transfe r, the method was 
first used by Mihala setall dl978l) . and later by dOlson & Kunaszl 



U987tlKunasz & Auerlll988HKunasz & oTsonl fT988). and subse- 
quently in many studies - see references below. 

3.1. Overview of the Method 

Let us consider, at time t, the radiation propagating in direction 
n. The optical depth from position r to position r + sn, where s 
is the path length between these two positions, is: 

r(r, r + sn, v, f) = I %(r + s'n, n, v, t) ds' (4) 
Jo 

For a radiation propagating from an upwind position r u to a 
current position r c , the integral form of the formal solution of 
the time-independent RTE is: 

I c = I ue - T "' + f K S(T)e- (T -- T) dr , (5) 
Jo 

where I C =I (r c , n, v, t), /«=/ (r„, n, v, t), n is the unit vector along 
the straight line (r„, r c ), t is the optical depth from r„ to an in- 
termediate position r„ + sn, with < s < \\r c - r„\\, and t uc is 
the optical depth from r„ to r c . 

We consider a three-dimensional Cartesian grid (see 
Figure [TJ defined in the (O, x, y, z) coordinate system with the 
unit vectors (e x , e y , e z ). The grid cells are rectangular boxes with 
irregular spacing in each direction. The sizes of the cells in 
each direction depend on the position in this direction: Ax c& \\(x), 
AyoeflOO, and Az ce u(z). The state parameters, i.e., the mass den- 
sity, the temperature, and the velocity components in the ob- 
server's frame, are defined at each grid point. They are pro- 
vided by the (radiation)(magneto)hydrodynamics (RMHD) sim- 
ulations. A direction of propagation of the radiation, defined by 
the unit vector n, is specified by the polar angle 6 between n and 
the unit vector e z , and by the azimuthal angle <p, between e x and 
the projection of n on the x-y plane. In order to span the whole 
4n sr angular domain, we impose 9 e [0, n] and <p e [0, 2jt[. 

The task is to calculate the specific intensity I c in a cur- 
rent point M c in direction n. M c is defined by the intersection 
of the cells marked with indices i, i+l in ^-direction, j, j+l in 
y-direction, and k, k+1 in z-direction (see Fig. [T). This is ac- 
complished by using Equation ©. Following the propagation of 
the ray that goes through M c in direction n, we define a short- 
characteristic by the line joining the intersection of the ray with 
the upwind cell, M u , to the current point M c (the upwind cell is 
defined by the indices (i, j, k) in Fig. [TJ. I c can be determined if 
we know the following quantities: the upwind specific intensity 
/„, the source function Sir, n, v, t) and the absorption coefficient 
X(r, n, v, t) (in order to deduct the optical depths t and t„ c ) along 
the short-characteristic from M u to M c . However, S and^f are 
specified only in the vertices of the cells (grid points). Therefore, 
we are essentially free to define laws of variation of these quan- 
tities along the short-characteristic, typically as low-order poly- 
nomials; we choose third degree polynomials. To do so, we need 
to know (S,x) in the upwind end point M u , (S u ,x u ), and in the 
downwind end point Mj, (Sd,X<t)- is the intersection of the 
short-characteristic with the downwind cell, which is defined by 
the indices (i + 1 , j + 1, k + 1) in Fig. [TJ We also choose third de- 
gree polynomials to get the upwind quantities (/„, S u ,Xu) and the 
downwind quantities (S d,Xd), through interpolations from the 
values in the neighboring grid points. Further, in Section [4] we 
describe and discuss in detail the discretization of Equation d5J, 
and the mathematical form of the third degree polynomials. 
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Fig. 1. Short-characteristics method illustrated with an example 
of a 3D irregularly spaced Cartesian grid. We consider all quan- 
tities in the observer's frame. They are defined in the vertices of 
the cells (grid points). (O, x,y, z) is an example of a coordinate 
system with the unit vectors (e x ,e y ,e z ). The specific intensity is 
calculated in the current point M c , for a radiation propagating 
from the upwind end point M u to the downwind end point Mj. 
The direction n of the ray is defined by the polar angle 9 and the 
azimuthal angle (p. The short-characteristic is defined by the line 
joining M u and M c . The transfer equation is solved in its integral 
form along this short-characteristic. The cells around point M c 
are marked with the following indices: i, i+\ in x-direction, j, 
7+1 in y-direction, and k, k+l in z-direction. 



3.2. Intersections of a Short-Characteristic with the 
Neighboring Cell Faces 

The coordinates of the intersections of a short-characteristic with 
the neighboring cell faces, i.e., the coordinates of the upwind 
end point M„ and the downwind end point Mj, are easily deter- 
mined in the local coordinate system defined by (M c ,e x ,e y ,e z ) 
(see FigJTJ- 

Figure [2] displays the upwind cell (top panel) and the down- 
wind cell (bottom panel), as they are defined for the short- 
characteristic plotted in Fig.Q] In this example, the upwind cell is 
determined by indices (i, j, k), while the downwind cell is deter- 
mined by indices (z+1, j+l, k+l). Let us refer to f x (respectively 
f y , f z ) as the generic name of a cell face that is perpendicular to 
x-axis (respectively y-axis, z-axis), and that may be intersected 
by a short-characteristic, i.e., to which M u or Mj may belong. 
In our example, f x in the upwind cell is specified by x — -Ax,-, 
while f x in the downwind cell is specified by x = Ax,+ i, both in 
the local coordinate system. In the same vein, f. is y — -Ayj in 
the upwind cell, and y - Ayj+i in the downwind cell. Finally, f 
is z — -Azk in the upwind cell, and z - Azk+i in the downwind 
cell. 

In the example displayed in Fig. Q]and[2] M u belongs to the 
face fz'. Z = -Azk, and Mj belongs to the face f z :z = Azk+\- 
However, depending on the values of (9, <p), M„ may belong to 
fy'.y- -Ayj, or to f x :x= -Ax,-, and Mj may belong to f y : 
y = Ay ;+ i, or to f x : x = Ax,+i. In addition, the upwind cell and 
the downwind cell are respectively (z, j, k) and (z + 1, j+l, k+l) 
only for (9,tp) e [0,7r/2[x[0, n/2[, as exemplified in Fig. Q] and 
|2] In fact, for each grid point, eight different upwind cells are 
possible, one different cell for (9, <p) in each of the eight octants 
in the An sr directional domain. Let us refer to Ax„ and Ax^ as 



face/, : z = Az k+S 



k + l 







Z ' 

























facef y : 
y = Ay J+ i 



face/ t 
x = Ax.u 



Fig. 2. Upwind cell (top) and downwind cell (bottom), for the 
short-characteristic M U M C displayed in Fig. \T\ In this example, 
the direction of propagation of the radiation, n, is in the first 
octant, which corresponds to (6, <p) e [0,7r/2[x[0, tt/2[. The 
local coordinate system is defined by (M c z ). The up- 

wind end point M u belongs to the cell face f z : z — -Azk- 
However, depending on (9, <p) in this octant, M u may belong to 
f x ■ : x — -Ax,-, or to fy'.y- -Ayj. In the general case, Ax, + Ayj, 
Ax, + Azk, Ayj ± Azk, though equalities are possible. The down- 
wind end point belongs to f z :z = Azk+i- However, depend- 
ing on (9, (p) in this octant, M c { may belong to f x :x = Ax, + i or 
to f . : y — Ay j+i. See Section [3~!2l for more details. 



the generic names for the size in x-direction of, respectively, the 
upwind and the downwind cell, for a given current point M c and 
a given direction (9, tp). Similarly, let us introduce the analogous 
notations Ay„, Ay^ and Az u , Azd for the corresponding sizes in y 
and z directions. Table [^indicates, for each of the eight octants, 
the corresponding range of (9, tp), the indices of the upwind cell, 
along with the values of Ax M , Ay M , Az u , for any current grid point 
M c . Table|2]indicates the equivalent quantities for the downwind 
cell. 
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The following generic inequalities identify the cell face to 
which an upwind (respectively downwind) end point belongs: 



M u(d ) e fx « |tan^| < ,|tan6»cos^| > 



M u(d ) £ f y <=> | tan > , | tan 8 sin tp\ > 



Ax 



u(d) 



^Zu(d) 



(6a) 
(6b) 



A^uCO € fz I tan cos ^| < , |tan0siny>| < — — (6c) 



The local coordinates of the upwind (respectively down- 
wind) end point M„(,/) depend on the face it belongs to, as fol- 
lows: 



M u{d) 



x U (d) - a xu (d)Ax U (d) 
y U (d) = u xu(d )Ax U (d) tan <p 

Zu(d) - a xu (d) tan g coslf 



M m : 

6/, 



v - ™ A - V " (rf) 

y«(d) = a yu (d)Ay U (d) 

_ Av„( rf ) 



= Az„ w tan cos (/> 
= a z «(d)Az MM tan0sin^ 

Z«(<i) = O zu (d)Az u (d), 



(7a) 



(7b) 



(7c) 



where the factors (a xu (d), a zu (d)), are equal to +1, depend- 
ing on the octant to which the direction n belongs, as detailed in 
Tables Sand 



4. Piecewise Cubic, Locally Monotonic, 
Interpolation 

4. 1 . The Advantages of a Piecewise, Locally Monotonic, 
Interpolant 



As mentioned in Section I3.ll solving the integral form of the 
radiative transfer equation with the short-characteristics method 
requires to define laws of variation of physical quantities within 
faces of the grid and along the short-characteristics. The lin- 
ear law is the simplest and fastest method. However, it leads 
to large numerical diffusion. This means that a sharp beam is 
significantly dispersed as it propagates th roughout a grid, which 
was discussed in detail by several auth ors ((Kunasz & Auer 1988; 
lHavek et al.l 120101 : iDavis et all l2012h . Moreover, a second or 
higher order polynomial is mandatory, in order to recover the 
diffusion approximation at large optical depth. A parabolic law, 
however, generates overshoots that may be negat ive and thus un- 
physic al. This point was discusse d in detail by iKunasz & Auerl 
( 1988]) and lAuer & Paletoul ( 1 1994b . An efficient method to over- 
come this difficulty is to use in each direction (x, y, z and the 
short-characteristic's direction) a piecewise conti nuous, locall y 
monotonic, interpolation function, as proposed bv lAuerl (|2003). 
The monotonicity suppresses any spurious extrema. Hereafter, 
we describe the method that we have adopted. 



MY 















W i 1/ 













Fig. 3. Piecewise interpolation. w(x) is the function to be in- 
terpolated between the end points (or nodes) jc, and Xi+\. We 
assume that we know w and its derivatives at the end points: 
Wi = w(Xi), w M = w(x,+i), w' t = w'(xi), w| +1 = w'(x,+i). The 
local interpolation function, H((x), is a cubic Hermite polyno- 
mial. It is a third degree polynomial defined between Xi and 
Xt+i, from the four relations: H(xi) = Wt, H(x, + i) = w,+i, 
H' (xi) - w'j, H' (xi + \) - w' i+v The arrows indicate the slopes 
of H(x) at the end points, or, equivalently, the derivatives w\, 
w'. 



See Section I4~2l for detailed explanations. 



positions or nodes (x,), = i „, with < X\ < x^ < ... < x n . In 
other words, w,- = w(xi) for i = 1,2, ...,n are given. Our objec- 
tive is to find an interpolation function that goes through all the 
data points (xt, wi) and that preserves the monotonicity of the 
data. A solution to this problem is provided by a piecewise cubic 
Hermite polynomial coupled with a weighted harmonic mean 
defining the derivative at each node, w' i = w'(xi). We detail this 
solution hereafter. 

Among the nodes above, let us choose two consecutive ones, 
Xi and Xt+i (see Figure [3). Within the interval [xt, Xt+i], a cu- 
bic H ermite polynomial is defined as (see equation (2) of I Auerl 
2003): 



H i (x) = {\-3q 1 i +2qf)w i 
+ (3q 2 i -2qf)w M 
+ [q] - 2q\ + q,j (x M - xi) w[ 



(8) 



where 



q t (x) 



x i+\ x i 



< qi(x) < 1 



(9) 



By definition, the Hermite polynomial matches the interpo- 
lated function w(x) and its derivatives at both ends of the inter- 
val: 

f Hi (xi) = Wi 
\Hi(x i+ i) = w M 



(10a) 



4.2. Cubic Hermite Polynomial and Weighted Harmonic 
Mean Node Derivative 

Let us consider a physical quantity that depends on only one 
variable, w(x), and which is specified only in a set of m discrete 



H' i (x i ) = w' i 
H'^Xj+i) = w\_ 



(10b) 



The derivative at an inner node, (w')j = 2. n -i, is define d by 
the weighted harmonic mean suggested by iBrodlid dl980h and 
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Table 1. Correspondence between octant, directional angles (0, (p), upwind cell indices, sizes of the upwind cell faces in each 
direction, and (a xu , a yu , a zu ) factors. (See Sections [3.1l and l3~2l for the definitions of these quantities.) 



Octant 


(0, <p) (radians) 


Upwind cell 
indices 


Face f x 
Ax 


Face 
Av 


Face f z 
Az 


ttza) 


1 


[0,f[x[Q,f[ 




1 - Ax ( | 


1 - Ay,| 


1 - Aa| 


(-1,-1,-1) 


2 


[0,f[x[f,7r[ 


(! + l,;,fc) 


Ax f+ i 


1 - Ay,| 


1 - A Z{ .| 


(+1,-1,-1) 


3 






Ax-, i 




1 - 


C+l +1 -1) 


4 


[0,f[x[f ,2;r[ 




1 - Ax,| 


A >V+i 


1 - Aa| 


(-i,+i,-i) 


5 


[f,*[x[0,f[ 




1 - Ax,| 


1 - Ay,| 


Azt+i 


(-1,-1, +i) 


6 


[f,7r[x[§,;r[ 


(i+ 1) 


Ax j+1 


1 - A>vl 


Azt+i 


(+1,-1, +i) 


7 


[f,>r[x[;r,& [ 


+ 1,&+ 1) 


Ax i+1 


Ay, + i 


Az^+i 


(+i,+i,+i) 


8 


[f.aM* 2»r[ 


(i,j+l,Jfc+l) 


1 - Ax,-| 


Ay;+i 


Azt+i 





Notes. We intentionally use the absolute value notations above, in order to emphasize the coordinates, in the local coordinate system, of the upwind 
cell faces that may be intersected by a short-characteristic. For example, for octant 1, the x-coordinate of f x is x = -Ax,. 

Table 2. Same as Table[TJ but for the downwind cell faces. 



Octant 


(6, ip) (radians) 


Downwind cell 
indices 


Face f x 
Ax rf 


Face f y 
Ay rf 


Face / z 
Az,i 


(a zd ,a xc i,a yd ) 


1 


[0.f[x[0,f[ 


i,fe + i) 


Axj+i 


A.V;+I 


Azt+i 


(+1,+1,+1) 


2 


[0,§[x[f,/r[ 


(U+U+l) 


1 - Ax;| 


Ay ; +i 


Az*+i 


(-1,+1,+D 


3 


[0,f[x[w,$[ 




1 - Axd 


1 - Ay y | 


Az^+i 


(-1,-1, +1) 


4 


[0,f[x[f ,2w[ 


(i+ 1, j,/t+ 1) 


Axj+i 


1 - Ay y | 


Az*+1 


(+1,-1, +i) 


5 


[f,w[xtaf[ 




Axi+i 


A>V+i 


1 - Az*| 


C+l, +-1,-1) 


6 


[f.w[x[f,w[ 


l,k) 


1 - Axd 


A.V;+1 


1 - Az,| 


(-1.+1.-1) 


7 


[f,w[x[jr,^[ 




1 - Ax,-| 


1 - Ay y | 


1 - Aza-I 


(-1,-1,-1) 


8 


[f,/r[x[f ,2w[ 


(i+hj,k) 


Axj+i 


1 - Ay y | 


1 - Az*| 


(+1,-1,-1) 



for x e]xj, Xj+i[, is equal to the sign of /.'. Therefore, spurious 
extrema never occur, which guarantees the positivity of the in- 
terpolated physical quantities. Let us note H(x) the piecewise 
interpolation function defined over the whole interval [x\, x„] = 
\J"=i[xi,Xi+i], so that its restriction to each interval [xy, xy+i] is 
equal the local function Hi(x). H(x) is continuous. Moreover, it 
is easy to demonstrate that its derivative H'{x) is also continuous 
throughout the whole interval [x\, x„], specifically at each inner 
node (x,),=2, ; ,-i- 

In addition to the Hermite interpolation. lAueil |2003) sug- 
gests a second possibility: using a Bezier polynomial as a lo- 
cal interpolant. The latter is quite close to a Hermite polyno- 
mial. Both functions match the local end point values w,, w,+i. 
However, while the Hermite polynomial matches the end point 
derivatives w' v the Bezier polynomial does not do it nec- 
essarily. The reason is as follows. Let us note B,(x) the local 
Bezier polynomial defined in [x;, x,+i]. fi,(x), and, therefore, its 
deriv ative ffl(x), depend on free parameters, called "control val- 
ues ' (lAuerll2003l) . These values can be adjusted in order to con- 
trol the variation of fi,(x) in the vicinity of the end points. In 
particular, they can be restrained within a minimum value and 
a maximum value (see equation (11) of lAuerl 12003). in order 



iFritsch & Butlandl d 19841), and brough t to the attention of the as- 
trophysics community by lAuerl (120031) : 



/j-i/j 



o if/'-,// so. 



where 



n = 



Wj+l - Wj 
x i+l — x i 



1 / X i+i - Xi 

a t = -11 + 

3 \ x M -Xi-\ 



(11a) 



(lib) 



(He) 



The derivatives at the outer nodes, w\ and w' n , are defined by 
the local slopes: 

(12a) 



x 2 - x\ 
w„ - W„-\ 
%n ~ %n— 1 



(12b) 



Adopting this definition of the node derivatives, each local 
cubic Hermite interpolant (H,(x)), = i. ;l _i is monotonic in the in- 
terpolation interval [x,, x, + i]. More precisely, the sign of Hj(x), 



6 



Ibgui et al.: IRIS, a generic 3D radiative transfer code 



to guarantee the monotonicity of Bi(x). In this case, the conti- 
nuity of the derivative B'(x) of the piecewise interpolant B(x), 
which is defined over [xi,x„] similarly to H(x), is not neces- 
sarily verified at the nodes. In short, whatever the definition of 
the node derivatives of w(x), (w')i=i,«> a piecewise Hermite in- 
terpolation function, H(x), guarantees that H'(x) matches these 
node derivatives and that H'{x) is continuous at the nodes. On the 
other hand, a piecewise Bezier interpolation function B(x), along 
with an adequate choice of the control values, guarantees to be 
locally monotonic, but B'(x) does not necessarily match the node 
derivatives and is not necessarily continuous at the nodes. Now, 
if we adopt very specific control values (see equation (9) of lAuerl 
2003 in the cubic case), we can force the Bezier polynomial to 
be identical to the Hermite polynomial. And, as we discussed in 
the paragraph above, with an adequate choice of the definition of 
the nodes derivatives (wj), = i in (Equations[TT1and[T2l. a piecewise 
Hermite interpolation function can be made locally monotonic. 

Finally, let us mention another possible definition of the node 
derivative s wU that ens ures the monotonicity of a Hermite poly- 
nomial. Steffen] (1 1990b suggests to calculate the slope at inner 
node (x,) I= 2,«-i from the unique parabola passing through the 
points (jtj-i, W;_r, Xi, Wi, Xt+i, Wj+i), if this parabola is monotonic 
in [x,_i, If not, the node derivative is defined as the com- 
mon slope shared by two parabolas, which are locally monotonic 
in [Xi~i, and [x,-,x !+ i]. Specific definitions are proposed for 
the outer nod es derivatives w[ and w' n . Although this possibility 
sugg ested bv lSteffeddl990l) satisfies our requirements of mono- 
tonicity, we did not implement it in IRIS, because the calcula- 
tion of such a derivative, for a given node, uses a larger cpu time 
than the calculation with the weighted harmonic mean formula 
(Equation (II lab in this paper may be compared to equation (11) 
of lStefFenlfl990i) . 

5. Implementation 

5.1. Interpolations in Cell Faces, Ghost Nodes 

For each direction of propagation of the radiation, we apply 
the 3D short-characteristics method, which consists in solv- 
ing the integral form of the RTE, Equation (0, throughout the 
3D Cartesian grid points. We repeat this procedure for all the di- 
rections specified over a solid angle. The method that we employ 
for propagating the radiation gives a privileged role to z-axis, in 
the sense that we are sweeping the 3D grid gradually, z-plane 
by z-plane, in the direction of increasing (or decreasing) z. In an 
astrophysical context, z-axis should represent a global direction 
of the energy transport in an object; that is the direction from 
the interior to the outer layers of a stellar or a planetary atmo- 
sphere, the vertical distance from an accretion disk plane, or the 
direction of an accretion column in a young stellar object. 

Boundary conditions must be defined in appropriate planes, 
which depend on the direction of propagation of the radiation. 
Assuming that the grid is defined by N x cells in x-direction, N y 
cells in y-direction, and N z cells in z-direction, the positions of 
the grid nodes in each direction may be noted as (xq, x\, x^ t ), 
(yo,yi, ... ,yN y ), and (zo,Zi,—,Zn z )- The sizes of the cells in 
each direction are then noted (Axi, Ax^J, (Ay\, Ay^ ), and 
(Azi, AzaO, and me y are defined by: 

<&yj = yj-yj-i (13) 
AZ* = Zk-yk-i 

In agreement with our notations in Fig. [TJ a propagation 
in the first octant corresponds to directional angles (0, <p) e 



[0, 7r/2[x[0, tt/2[, and, therefore, to a propagation along increas- 
ing x, increasing y, and increasing z. Consequently, the bound- 
ary specific intensities must be known at the bottom frontiers of 
the 3D Cartesian grid, defined by the following planes x = xq, 
y = yo,z = zo- 

As explained in Section |3~T1 calculating the formal solution 
of the monochromatic RTE in its integral form, at a given grid 
point and for a given direction of propagation of the radiation, 
requires the determination of (I u ,S u ,Xu) at the intersection of 
the short-characteristic with the upwind face, and {S d,Xd) at the 
intersection of the short-characteristic with the downwind face. 

The opacities x are known at each point of the three- 
dimensional grid. And, by definition of a formal solution, the 
source function S is known at each grid point. Then, S„, Xm 
(respectively S4, xi) are determined by successive cubic mono- 
tonic interpolations in each of the two directions of a given up- 
wind (respectively downwind) horizontal of vertical face. The 
top panel of Figure |4] provides an illustration of a sequence of 
interpolations in a horizontal upwind cell face. This example 
corresponds to the case shown in the top panel of Fig. [2] The up- 
wind face is identified by indices i and j and belongs to the plane 
defined by f z : z — -Azk in the local coordinate system (see 
Section l3~2l . The ray propagates in the first octant, i.e., along in- 
creasing x, y, and z. The black dots indicate the grid points. M„ 
is the upwind end point. M^r is the projection of the current point 
M c on f z , and coincides with a grid point. Quantities defined in 
the grid points are interpolated in y-direction at the y u position 
of M u . Then, the interpolated values, in positions A;_2, Aj_i, A,, 
Ai + i, are themselves interpolated in x-direction at the x„ position 
of M u , which provides the upwind quantity S u , Xu- The curved 
arrows indicate the grid points which contribute to the interpola- 
tions. Four points are necessary to accomplish each interpolation 
(see Section l4~2l along with Equations <[8J and 03} )- 

Unlike opacities and source functions, the upwind specific 
intensity is not always known in all the grid points neighbor- 
ing the upwind cell face. The sweeping procedure throughout 
the 3D Cartesian grid is made z-plane by z-plane, as described 
above in the first paragraph of this section. Therefore, for each 
new z-plane to be processed, the specific intensity is known at 
each grid point of the upwind z-plane. Accordingly, the proce- 
dure described above for S andx works if M u belongs to a hori- 
zontal face. However, the specific intensity is not known at all 
grid points of a preceding vertical plane (j x : x = ±Ax,-, or 
f y : y = +Ayj, in the local coordinate system), except if it is 
a boundary plane. The bottom panel of Fig. [4] provides an ex- 
ample of a sequence of interpolations of the specific intensity 
in a vertical upwind cell face. In this case, z = Zk is the current 
z-plane, z = Zk-i is the upwind one. M„ belongs to the vertical 
face identified by indices j and k, and by f x :x- -Ax,. We still 
consider a ray propagating in the first octant (increasing x, y, z). 
The specific intensity is fully known in the preceding z = Zk-i 
and z = Zk-\ planes, where we can perform interpolations in y- 
direction at the y u position of M u . This provides the values of the 
specific intensity in points A^-2 and At_i. The crosses indicate 
the grid points in which the specific intensity is not known. Since 
the sweeping follows the direction of propagation of the ray, it is 
performed along increasing y, in the z = Zk plane. Therefore, the 
specific intensity is known in the three consecutive grid points 
((yj-2,Zk), (yj-uZk), (yj,Zk)), but not in (yj+uZk)- 11 is also not 
known in the z = Zk+i plane (except in y — y ; -_2 if this y-plane 
is a boundary). Consequently, we can determine the specific in- 
tensity in point A^ in z = Zk plane, with a cubic Hermite inter- 
polation in y-direction, using the values in ((yj-2, Zk), (yj-i,Zk), 
(yj,Zk))- Finally, the interpolated values in At_2, An, Ak are 
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themselves interpolated in z-direction at the z u position of M u , 
which provides the upwind specific intensity /„. For the last two 
interpolations, since we know the intensities in only three grid 
points, the derivatives of the cubic Hermite polynomials in the 
edge points, M^r in y-direction and Ak in z-direction, are pro- 
vided by the local slopes as defined in Equations (fTSt . 

We have added at each boundary of each of the three x, y, 
z dimensions, two layers of ghost nodes that embed the initial 
3D Cartesian grid. These ghost nodes are necessary if we want to 
apply the above interpolation procedure uniformly to all the grid 
points, which spares us from testing whether the grid points are 
close to the boundaries or not. In ^-direction, for example, and 
for a propagation in increasing x, ghost nodes are required for 
determining the upwind quantities at X\, and for determining the 
downwind quantities at jc^-i and xn . The ghost nodes (X-%, 
and (Xff M , Xff M ) are added. Their positions per se do not matter, 
as we define in these nodes positive linear extrapolations of the 
state parameters. Then, it is easy to show that we can still use 
the cubic Hermite polynomial interpolants in x\, Xj^ x -\ and xn x , 
and that the derivatives in these nodes, even though calculated 
with the weighted harmonic mean, are equal to the local slopes 
defined by Equations (fT2l >. 



5.2. The Discretized Integrated RTE along a 
Short-Characteristic 

In the preceding subsection, we have described our procedure to 
determine the upwind quantities (I U ,S u ,Xu), and the downwind 
quantities (Sd,Xd)- Once we know them, we are in a position to 
calculate the integral form of the RTE, Equation (0. 

To do so, we define laws of variation, along the short- 
characteristics, of the source function and the opacity. We use 
a monotonic cubic Hermite polynomial between the upwind 
endpoint M u and the current grid point M c (see Fig. [TJ. Since 
we know the values of S and \ m three points, M u , M c , and 
the downwind endpoint Mj, the derivatives of these quantities 



at M„, ^r| H and ^| , are provided by the slopes defined by 
Equation ( 112ab . while the derivatives at M c 



ds_\ 

dr I 



and 



m 

ds 



are 



calculated with the weighted harmonic mean, Equation ( fTTT ). 
Then, the RTE (Equation (|5]l) can be integrated analytically. We 
obtain the following expression: 



I c = I u e T,IC +aS u +pS c + a' 



dS 
(It 



dS 
dr 



where 



a = — {6 (juc - 2) + [-4 + 6 (t uc + 2)]e~ T "> 



P = —{[rl - 6 {t uc - 2)] - 6 (t uc + 2) 
a' = 4-{(2t„ c - 6) + (tI + 4r„ c + 6) 

' 'uc 

P = -^{{-rl + 4r„ e - 6) + (2r„ e + 6) 



(14a) 

(14b) 

(14c) 
(14d) 
(14e) 



The optical depth t uc is calculated by integrating analytically 
Equation © from M u to M c . We obtain: 



- 1 1 2 d X 

T uc — ^ S " c OCu + Xc) + 2 S uc I 



u ds 




z k+\ X — X-- 




y } -2 yj-i Ak - 



'y 



(15) 



Fig. 4. Examples of cubic monotonic interpolations in cell faces 
intersected by a short-characteristic, for a ray propagating in the 
first octant, i.e., along increasing x, y, and z. Two configurations 
are possible. The first one is displayed in the top panel, where all 
the quantities are known in the grid points neighboring the cell 
face. This situation always occurs for the upwind and downwind 
source functions and opacities (S u ,x u ), (Sd,Xd)- But, for the up- 
wind specific intensity /„, it occurs only when horizontal faces 
f z : z = +Azk or boundary faces are intersected. The second con- 
figuration is displayed in the bottom panel. It impacts the deter- 
mination of /„ when vertical non-boundary faces, f x : x = +Ax,- 
or f y : y = +Ayj, are intersected. In this case, the upwind spe- 
cific intensity is not known in all the grid points neighboring the 
cell face. See Section l5Tl for detailed explanations. 



where s uc is the path length from M u to M c . 

When Tuc is small enough (< 5 x 10~ 2 ), we use Taylor expan- 
sions up to the third order of Equations (fT4l . The bottom line is 
that the RTE can always be calculated, including the two asymp- 
totic limits: the optically thin or even transparent path for which 
Tuc equals zero, and, therefore, I c — > /„, and the optically thick 
or even opaque path for which t uc tends to infinity, and, there- 
fore, — > S r - 



5.3. The Velocity Gradient Effect (Doppler shift) 

We consider a medium with a non-zero macroscopic velocity 
with respect to the observer's frame. In this case, the velocity 
gradient between two positions in the medium causes a Doppler 
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shift of photon frequencies, and therefore a shift of the cor- 
responding spectral lines. This effect must be properly taken 
into account in an implementation of the short-characteristics 
method. If the velocity difference between the current point M c 
and the upwind end point M u (see Fig.[TJ is large enough, then a 
simple interpolation scheme may underestimate or overestimate 
opacities. 

To illustrate this point, let us consider a single emission line 
and a given frequency at which we solve the RTE. The following 
situation may happen. At M u , the frequency is located at a wing 
of the line, quite far from the center of the line, such that the 
absorption coefficient is very small compared to the line center 
value. And, at M c , the line is shifted, such that this same fre- 
quency is at the opposite side of the line, in the other wing, far 
from the center. The law of variation then infers a very weak 
opacity for this line throughout the path from M u to M c . In fact, 
the line goes through its maximum value at this frequency, some- 
where between M u and M c , a nd this is not taken into account. 
This issue was pointed out by Ivan Noort et al. who em- 

phasized the fact that the opacity along a short-characteristic 
may be underestimated by several orders of magnitude. In ad- 
dition to that example, we provide here a case of overestimation. 
Let us consider now two lines that are quite distant from each 
other. We may face the following situation. At M u , the calcu- 
lation frequency can be located close to the center of one line. 
And, at M c , due to the shift of the lines, this frequency can be 
close to the maximum of the second line. Consequently, a vari- 
ation law defined between M u and M c overestimates the opacity 
if the two line centers are distant enough, such that we miss the 
trough between them. 

A commo n solution of these pro blems is provided by a sub- 
griding (e.g. . Ivan Noort et a Let v be the macroscopic 
velocity in the observer's frame, and v = v • n its projection on 
the direction of photon propagation n. The idea of subgriding is 
to subdivide a short-characteristic into a set of subintervals, so 
that the difference in v between two points of this subdivision 
remains small compared to the thermal velocity. In IRIS, we fol- 
low this procedure. Specifically, we introduce a dimensionless 
parameter factor e„ that provides to the user a control over the 
subdivision. Let us identify with indices I and / + 1 two points 
of the subdivision of a given short-characteristic. Let us con- 
sider a single line. We note vo the transition frequency (emission 
or absorption of a photon) in the particle's frame (or comoving 
frame), vo coincides with the center of the line in the fluid rest 
frame where v = 0, but this frequency is shifted to the value v ct r 
in the observer's frame where v + 0, as follows, 



Vctr = V 



(16) 



where we consider a Doppler shift of frequencies to first order 
in v/c, since we assume the velocities to be non relativistic (see 
Section |2]i. v ctl is the center of the line in the observer's frame. 
The shift of v ctr between position / and position / + 1 is given by 



v c tr(f + l)-v ctr (/) = v 



v(/+l)-v(/) 



(17) 



The Doppler width associated with this line, in a given position 
in the medium, is 

Av D = v — (18) 

c 

with c being the speed of light, and V t h the thermal velocity. The 
latter is defined as 



Vth 



2kT 



1/2 



(19) 



where k is the Boltzmann constant, m the mass of the particle. e D 
is a control parameter defined such that the Doppler shift of the 
center of the line between position / and position / + 1 is bounded 
as follows: 



r(/+l) 



x(OI < e D 



Av D (0 + Av B (Z+l) 



(20) 



The quantity j[Av D ([) + Av D (l + 1)] is an average of the local 
Doppler width between its value in I position and its value in 
/ + 1 position. It is easy to demonstrate that this inequality is 
independent of the transition frequency in the particle's frame 
vo, and that it is equivalent to the following condition on the 
gradients of the projected velocities: 



|v(/ + 1) - v(OI < e„ 



Vth(/) + Vth(/+1) 



(21) 



e D is a tunable parameter, whose value depends on the required 
accuracy of the solution (ideally between 1/3 and 1). The smaller 
it is, the higher accuracy is achieved; however at the expense of 
increased computer demands. 

6. Radiation Moments and Angular Integration 

In addition to computing specific intensity, one can calculate im- 
portant angle-averaged quantities related to radiation moments, 
such as the mean intensity J(r, v, t), the radiation flux vector 
F(r, v, f), and the radiation pressure tensor P(r, v, t). These are 
defined as 

J(r, v, t) = — £ I(r, n, v, t) dQ, (22) 



An 



F(r, v, t) 



(£) I(r, n, v, t) n 



dQ, 



P(r, v, t) = - I(r, n, v, t) nn dQ., 



(23) 
(24) 



where dQ. is the elementary solid angle, and the symbol j> desig- 
nates an integration over An sr. The flux vector and the radiation 
pressure tensor can also be written in component form: 



Fj = (£) I(r, n, v, t) m dQ, 
P u = - ^) I( r , n, v, t) riitij dQ, 



(25) 



(26) 



In Cartesian coordinates, the three directions (1,2,3) are the unit 
vectors (e x , e y , e z ) introduced in Section l3~Tl and the components 
of the radiation propagation vector n are 



n\ — n x = n ■ e x = sin cos tp 
«2 = n y — n ■ e y = sin 6 sin tp 
«3 = n, = n ■ e z = cos 6 



(27) 



IRIS computes each of the three components of the radiation flux 
vect or, F x , F v , F 7 . Since the radiation pressur e tensor is symmet- 
ric dMihalaslll978l:lMihalas & Mihalasll 19841) . it is defined by six 
components, P xx , P vy , P zz , P„, P xz , P yz , all of them being com- 
puted by IRIS. 

The angular integrations are numerically performed by us- 
ing quadratures. A quadrature is a set of M discretized direction 
vectors n m = (n mx , n my , n mz ) and the corresponding weights w m , 
which are the numerical representation of, respectively, the di- 
rection vector n and the elementary solid angle dQ. An appro- 
priate quadrature for an integration over the whole 47r sr range 
should satisfy at least the first two of the three following condi- 
tions: 
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1. symmetry: the directions of the quadrature are invariant af- 
ter rotation of 90° around the three coordinate directions 



for r 



2. for an isotropic radiation field, the three moments should be 

numerically exact. 

3. for an isotropic radiation field, the half first moment should 

be numerically exact. 
The second condition reads as follows: 



M 



dD, -An - 2^ w m, 

m=\ 
M 

ndQ, = = ^w m n m , 
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An \— i 
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(28a) 
(28b) 
(28c) 



where 6 is the unit tensor. The third condition is optional, and 
is related to the following issue. One may need to calculate a 
radiative flux over a surface or at a wall. In this case, the flux is 
determined by an integration over 2n sr, instead of An sr: F wa u = 
/ n dCl. The third condition reads then: 



2k 



I n dQ, = — | n dQ. = net — ^ i 

Jne t >0 Jnc,<0 n , c > q 



w m n„ 



(29) 



where et stands for e x , e y , or e z , and where we consider surfaces 
or walls perpendicular to one of these three directions. 

We h ave implemented in IRIS the quadratures of type A pro- 
posed by ICarlson| d 1963b. and w hose construction principle is 
summarized by Bruls et al. ( 1999). They satisfy conditions 1 and 
2, but not co ndition 3. We also have i mplemented the quadra- 
tures built bv lLathrop & Carlsonl ( 1 1965b . The latter provide sev- 
eral sets of quadratures that satisfy conditions 1 and 2, but not 
condition 3, with 24 up to 288 directions. They also provide sev- 
eral sets of another type of quadratures that satisfy the three con- 
ditions above, with 24 up to 80 directions. 

We also hav e implemented the Gaussian quadrature 
dPress et al.lll992h . The latter is not appropriate for the calcu- 
lation of the radiation moments, because it does not satisfy any 
of the three conditions above and leads to serious non physical 
asymmetries, unless one chooses a large number of directions 
(typically, more than 100). However, it is suitable for the inte- 
gration of the incoming specific intensity over a viewing solid 
angle, in order to determine the radiative power per unit surface 
that receives a detector looking at a given object. 

7. Periodic Boundary Conditions 

We consider a three-dimensional medium with an infinite exten- 
sion in the horizontal plane (x,y), and a finite extension along 
the vertical z-axis between its lower boundary zo and its upper 
boundary zn-- We assume that this medium has a double peri- 
odicity, one in x-direction, and one in y-direction. The bound- 
ary conditions are known at zo and Zn z - For example, we may 
consider a non-irradiated stellar atmosphere with no incoming 
radiation at the outer surface z# z and a black body radiation at 
its inner surface zo- The computational grid ranges from zo to 
Zn z in vertical direction, from Xq to xjf , and from yo to in 
the horizontal plane, so that (x^ x - xq) defines one x-period and 
OX ~yo) defines one y-period. Now, the 3D short-characteristics 
method consists in solving the integral form of the RTE by prop- 
agating the rays throughout a computational domain, from up 



to three upwind boundary sides in which the specific intensity 
is assumed to be known, down up to the three other faces of 
the domain. For horizontally periodic media, while the vertical 
boundary conditions are specified explicitly, the lateral bound- 
ary conditions are defined implicitly, such that, for any physical 
quantity f(v, x, y, z), we have the following relations: 

f(v, x Nx ,y, z) = /(v, x Q , y, z), for any y, z, (30a) 
f(v, x, y Ny , z) = f(v, x, y , z), for any x, z. (30b) 

Consequently, when the upwind end point of a short- 
characteristic intersects a lateral boundary, we prolong this char- 
acteristic, which becomes a long-characteristi c, until it intersects 
a h orizontal face, following th e suggestions bv lAuer et al 1 (fl994h 
and lFabiani Bendichol d2003i) . For a given direction of propaga- 
tion, and for each z-layer, this treatment affects only the rows 
that are adjacent to the boundary faces. 

Figure [5] illustrates our point with an example showing a 
propagation of the radiation in the first octant, i.e., along in- 
creasing x, y, and z. Two characteristics are plotted from the 
current point M c in the plane x = X\, which is adjacent to 
the boundary plane x — xq. The upwind end point M u of the 
short-characteristic (M u , M c , Md) belongs to a horizontal face. 
Therefore, the specific intensity can be calculated by interpola- 
tions from the values at the neighboring grid points in the up- 
wind horizontal plane. If we consider the short-characteristic 
(M' b , M c , M' d ), the upwind end point M' b belongs to the horizontal 
boundary plane x = Xq, in which the specific intensity is not ex- 
plicitly defined. The short-characteristic is then lengthened, and 
becomes a long-characteristic, down to its first intersection, M' u , 
with a horizontal face. The long-characteristic can go through 
several cells (only one cell in the example, for the clarity of the 
figure) beyond the vertical boundary face, before hitting a hori- 
zontal face. 

When the medium is not periodic, we start the sweeping of 
a given z-plane at the grid point that is the closest to the vertical 
boundary planes. For example, if the ray propagates in the first 
octant, these planes are x = Xq and y = yo, therefore the first 
grid point in which we calculate the specific intensity is (x\,y{). 
Then, the specific intensity is calculated progressively in all the 
grid points in the z-plane along x-direction and y-direction. Now, 
as we consider in this section a medium with a double horizontal 
periodicity, we can start the sweeping at any grid point in the z- 
plane, provided that the grid points where the specific intensity 
is calculated span one x-period and one y-period. We mark this 
starting point as (x !o ,y JO ). The indices z'o and jo are defined so as 
to minimize the number of cells that are intersected by the long- 
characteristics, thus saving computing time. This way, ;'o is cho- 
sen such that the corresponding upwind cell has the largest size 
in x-direction, Ax, = Ax max - Similarly, we choose jo such that 
Ayj = A y max . Then, starting from (x, , yj ), we calculate the spe- 
cific intensity with the long-characteristics method at y ;o , for the 
grid points along x-direction that span one x-period. In the same 
vein, we use the long-characteristics method at x, , for the grid 
points along y-direction that span one y-period. Once these two 
first rows have been treated, we resume the short-characteristics 
method for all the next rows of the current z-layer. Note that 
here we do not introduce ghost nodes (see Section lBTTl l. since the 
medium has an infinite extension in x-direction and y-direction. 

The schematic diagram of Figure [6] clarifies our point, for a 
radiation propagating in the first octant, therefore along increas- 
ing x, y, and z. It shows two computational domains in a given 
z-plane. Both span one x-period and one y-period. The original 
domain is depicted by the rectangle with solid thick lines. The 
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positions of the grid nodes in the horizontal z-plane were intro- 
duced in Section lBTTl (xq, X\, XnJ and (yo,y\, • ■■,yA< v )- The ver- 
tical boundary planes are x — xq andy = yo. The indices ;'o and jo 
are defined above in the preceding paragraph. The RTE is solved 
in the computational domain depicted by the rectangle with thick 
dotted lines. The positions of the grid nodes of this domain are 
(xi ,x io+u ...,x h+N J and (y jo ,y jo+u ...,y jo+Nv ), with the following 
correspondence of the grid sizes, derived from the double peri- 
odicity of the medium: 

A*/,/e[[A' l + lJv i+ / ] = Ax,v,,v e |[ Uo ]] 

We also have the following correspondence for the physical 
quantities /(v, x,y,z): 

fvAXuyj) ieWt+hNx+kl ~ fv^\Xi', yAfsllM (32) 

We employ the long-characteristics method along x-direction 
from (x k ,y jo ) to (x io + Nx ,y jo ), and along y-direction from (x, ,y ;o ) 
to (xi ,yj 0+ N r )- Then, we resume our usual short-characteristics 
method for the next rows of the current z-layer. Once calcu- 
lated, the specific intensity is defined in the original computa- 
tional domain, with a simple rearrangement of indices as per 
Equation d32l . Note that some long-characteristics may con- 
sist merely of short-characteristics. This is the case for a given 
grid point (x;,feR>,Ayi>:>Vo)> if the corresponding Ax t and Ay max 
are large enough, so that the related upwind end point inter- 
sects the horizontal upwind face. Similarly, this is the case for 
(xi >yj.Mh.N r f), if Ay ; - and Ax max are large enough. 

We mention here that it is also possible to handle peri- 
odic conditions by ite rativ e methods, as done for instance by 
Ivan Noort etall d2002h and iDavis et alJ d2012l) . However, such 
an approach can lead to a sl ow convergence i n the c ase of opti- 
cally thin media, as noted by Ivan Noort et a i1 (l200l . 

8. Tests of IRIS 

The code underwent a number of various tests to validate each of 
its functionalities that we describe in this paper. We focus here 
on the three most important ones: the searchlight beam test, a 
comparison with a well-tested one-dimensional scheme, and the 
test of the velocity gradient effect applied on one given spectral 
line. 

8.1. The Searchlight Beam Test 

In the context of astrophysical radiation transport, this test was 
first proposed by iKunasz & Auerl (I1988I) . Several authors then 
used it to eval uate their m ultidimensi onal short-chara ct eristics 
algorithms dStone et alJ Il992t lAuer & Paletoul [ l994; 
Fabiani Bendicho 2003; Fabiani Bendic ho & Truiillo Buend 
2007t lHavek et all l20ld IDavis et alJ I2012I) . The purpose is 
to examine how a beam propagates through a transparent 
medium. Theoretically, the beam crosses the medium without 
being absorbed or dispersed. The RTE, Equation (0, is greatly 
simplified since, in this particular case, I c — /„. Note that the 
numerical counterpart of the RTE, Equation (fT4l) is consistent 
with this behavior when the opacity is zero, as pointed out in 
Section I5.2I Therefore, the calculation uses only the upwind 
interpolation of the specific intensity. This test challenges 
the capabilities of the piecewise cubic, locally monotonic, 
interpolation technique (see Section @), as applied to cell face 
interpolations (see Section lBTTl ). 
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Fig. 5. Extract of a 3D Cartesian grid. Unlike the general case 
presented in Fig. [TJ we consider here the particular case of a 
medium with a double horizontal periodicity, in x-direction and 
in y-direction. In x-direction, one period is defined from xn to 
Xn x - The plane x — xq is a boundary plane. The cell defined by 
indices (—l,j,k), drawn with thick lines, is identical to the cell 
(N x , j, k) (not drawn). The upwind specific intensity is known at 
M u for the short-characteristic (M„, M c , Mj), since M u belongs 
to a horizontal face. On the other hand, the upwind specific inten- 
sity is not known at M' b for the short-characteristic (M' b , M c , M' d ), 
because M' b belongs to a vertical boundary face. Therefore, the 
characteristic is prolonged up to the first intersection with a hor- 
izontal face, M' u . See Section|7]for more explanations. 



For comparison, the numerical test was performed wi th the 
same parameters as the ones used by Hay ek et alJ (2010). We 
consider a hard-edged beam that propagates throughout a three- 
dimensional zero opacity box along a slanted direction defined 
by B = 28.1° and tp = 45.0°, as shown in the sketch (top left 
panel) of Figure [7] The box is made with 100 3 grid points. It 
extends in each of the three directions (x, y, z) from to 10 in ar- 
bitrary units. The beam is made with 30 2 grid points at the base 
of the box. It enters the box at the base and emerges from it at the 
top. The bottom panels of Fig.|7]show the normalized specific in- 
tensity of the beam as a function of x and y, at the base of the 
box (left panel), and at the top of the box (right panel). We ver- 
ify that the beam profile is very well conserved. There is no ab- 
sorption: the maximum value of the normalized specific intensity 
remains equal to one. The dispersion is very limited: the hard- 
edges are conserved, along with the size of the beam. In addition, 
the beam exits at the right position of the box. The positions of 
the four corners of the beam, [(xi,yi); (x2,yi); (xi,y2); (x2,y2)], 
are at the base: [(1 .5, 1 .5); (4.5, 1 .5); (1 .5, 4.5); (4.5, 4.5)], and at 
the top: [(5.3, 5.3); (8.3, 5.3); (5.3, 8.3); (8.3, 8.3)], which is con- 
sistent with the (9, tp) values specified above and the size of the 
box along z direction. In order to show in greater detail the dis- 
persion of the beam, we have plotted sectional views of the beam 
along x-axis (top right panel). The red curve represents the nor- 
malized specific intensity as a function of x, at z = (base of the 
box) and y = 3 (middle of the beam). The green curve represents 
this quantity at z = 10 (top of the box) and y = 6.8 (middle of 
the beam). In order to compare the two profiles, we have plotted 
the blue curve, which is the green one artificially shifted so that 
its center fits the center of the red curve. Such a superposition 
shows the symmetry of the beam's dispersion, which also re- 
mains small. We verified that our scheme guarantees the photon 
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Fig. 6. Original domain (solid thick lines) and calculation do- 
main (dotted thick lines) at a given z-plane, for a medium that is 
periodic in the horizontal (x,y) plane, and for a radiation propa- 
gating in the first octant, i.e., along increasing x, y, and z. Both 
domains span one period in x-direction and one period in y- 
direction. The nodes of the original domain are (xq,x\, ...,Xiy x ) 
and (yo,yi, —,yN y ). The nodes of the calculation domain, where 
the specific intensity is computed, are (Xi ,Xi +i,...,Xi +ff t ) and 
Cy/o'X/o+i' — >yjo+Ny)- The indices z'n and jo are chosen such that 
Axi„ = Axmax and Ay;„ = Ay max . See text in Section [7] for a dis- 



vi„ - Ax max and Ay jo 
cussion. 



conservation: the xy-surface integral of the specific intensity at 
z — equals the the xy-surface integral at z = 10, to the machine 
accuracy. 

This test represents an excellent validation of the piecewise 
cubic, locally monotonic, interpolation scheme that we use. Such 
a technique almost suppresses the numerical diffusion effect of 
the short-characteristics method, and is more efficient than lin- 
ear, parabolic (the first one is extremely diffu sive, the second 
one i ntroduces spurious extrema, as shown by i Kunasz & Auei 
1988), or even monotonic quadratic schemes dAuer & Paletou 
1994). 



8.2. Comparison with 1D Plane-Parallel Models 

It is important to check whether IRIS can reproduce the results 
provided by one-dimensional radiative transfer schemes, when 
the medium is made of homogeneous plane-parallel layers, i.e., 
in the case of a ID plane-parallel structure. Simulating the ra- 
diative transfer for such a medium with IRIS can be achieved by 
imposing periodic boundary conditions in faces perpendicular to 
the parallel layers. Assuming that these layers are parallel to the 
horizontal z-planes, the specific intensity has then the following 
dependence, I(z, cos 6, v). Applying the full 3D calculation with 
IRIS for such a medium perfectly reproduces the radiative re- 
sults provided by ID radiative models at any frequency, as we 
show in the example below. 

We calculate the radiative transfer through a radiative shock 
structure made of homogeneous plane-parallel layers. The hy- 
drodynamics structure, provided by Matthias Gonzale2Q, is ob- 
tained with th e three-dimensional r adiation hydrodynamics code 
HERACLES dGonzalez et alj2007h . The test case is a simulation 



of an experimental radiative shock generated in a tube full of 
Xenon assumed to be a perfect gas, with the following upstream 
conditions: fluid velocity = 60 km s _1 , pressure = 7 bar, temper- 
ature = 1 eV. The opacities that we use are from iMichaut et aD 
(120041) . Our objective here is focused on the comparison of the 
ID and the 3D radiative models; a future paper (Ibgui et al. 2012, 
in preparation) will present more details of three-dimensional 
models of radiative shocks. 

The 3D computational grid is built with 20 x 20 x 5 10 points 
in respectively x, y, and z directions. In each horizontal z-plane, 
all the state parameters (temperature, density, velocity) are inde- 
pendent of x and y. Two post-processing calculations were per- 
formed: one with IRIS applied on the 3D grid, and one with a 
well-tested ID short-characteristics solver applied on a ID grid 
that is extracted from the 3D one. Both were done at seven- 
teen frequencies that range from hv — 2 eV (A — 620 nm) to 
hv = 494 eV (A = 2.41 nm). The velocity gradient effect was 
not taken into account in this test. Figure [8] shows several radia- 
tive quantities provided by IRIS and by the ID solver. By way 
of example, we select one frequency, hv - 296 eV, which corre- 
sponds to A — 4.19 nm. 

The top left panel of Figure [8] displays the specific in- 
tensity as a function of z, and for different polar angles ft 
0° (green), 45° (blue), 89° (red), 135° (orange), 180° (cyan). 
The solid lines are the results calculated with IRIS. The 
symbols are the results provided by the ID solver. Note 
that, for the readability of the figure, we plot the ID re- 
sults for only some of the 510 z values. Note also that we 
plot the 3D results calculated with several azimuthal angles 
<p(°) = (0,45,89,90, 135, 179, 180,225,269,270,315,359). All 
the curves for these twelve <p angles are perfectly superposed, 
and, therefore, indistinguishable. This demonstrates that the spe- 
cific intensity is independent of <p. In addition, we verified that 
for any given pair (9, <p), the specific intensity at a given z-plane 
is the same for all the grid points of this plane, which demon- 
strates that this quantity is independent of x and y. Last but not 
least, the ID symbols perfectly fit with the 3D curves. 

The other three panels of Figure [8] display the moments ob- 
tained by angular integration (see Section|6]i: the mean intensity 
/ in the top right panel, the components of the radiation flux 
vector, F x , F y , F z , in the bottom left panel, and the components 
of the radiation pressure tensor P xx , P yy , P zz , P xy , P xz , P yz , in 
the bottom right panel. All the curves show that the ID and 3D 
results perfectly agree. We checked this result for all the seven- 
teen tested frequencies. In addition, the 3D results provided by 
IRIS verify, for any frequency, the following properties that are 
theoretically valid for ID plane-parallel structures. The radiative 
flux is zero in x and y direction, F x = 0, F y = 0, and it de- 
pends only on z in z-direction, F z (z, v). The non-diagonal com- 
ponents of the radiation pressure tensor are all zero: P xy = 0, 
Pxz = 0, P yz = 0. The other three pressure components depend 
only on z, P xx (z, v), P yy (z, v), P u {z, v), and we have, for any z and 
v,P xx {z,v) = P yy {z,v). 

These tests demonstrate that IRIS is capable of reproducing 
ID plane-parallel simulations, and, thereby, of handling periodic 
boundary conditions. 

8.3. Test of the Velocity Gradient Effect 



1 AIM, CEA/DSM/IRFU, CNRS, Universite Paris Diderot, 91191 
Gif-sur-Yvette, France 



As explained in Section 15.31 IRIS can handle the velocity gra- 
dient effect on spectra (caused by the Doppler shifts of pho- 
ton frequencies) by subdividing a short-characteristic in a set of 
subintervals. This subdivision is controlled by a tunable parame- 
ter e D , so that the Doppler shift of any spectral line between two 
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Fig. 7. Searchlight beam test. The top left panel of the figure depicts the purpose of this test. A slanted beam, with direction n, 
enters a zero opacity box at its base, and emerges from it at the top. The objective is to compare the entering beam profile with the 
emerging one; both are theoretically the same. The bottom panels show the normalized specific intensity of the beam, as a function 
of positions x and y, at the bottom of the box (left panel), and at the top of the box (right panel). The top right panel displays 
sectional views of the beam along x-axis, at z = and y — 3 (red curve), at z = 10 and y = 6.8 (green curve). The blue curve is 
identical to the the green one, but shifted so that its center fits the center of the red curve. See Section lOl for explanations. 



subinterval points remains bounded by a fraction e D of the local 
Doppler width of the line, as per Equation ( f20b . 

We illustrate the effect of our treatment in a simple ideal 
case that can be otherwise calculated quasi-analytically. We con- 
sider one plasma layer with a uniform temperature T in a stellar 
wind region with a velocity gradient, as sketched by Figure |9l 
top panel. We focus on one velocity direction, referred to as 
z-direction. In our example, the material moves toward an ob- 
server. The layer is numerically modeled with a single cell be- 
tween two grid points, z\, the farther position from the observer, 
and zzt the closer position to the observer. The wind velocities at 
grid points are V\ > and Vz > V\ . We assume that there is only 
one radiating species. We focus on the radiation of the layer it- 
self, in the direction of the wind toward the observer, and neglect 
any incoming radiation in z\- Although the situation in a stellar 
wind is significantly different from LTE, we do assume LTE for 
the purpose of this test. We also neglect scattering. Therefore, 
the source function is uniform and equals the Planck function at 
the layer's temperature, B v (T). We consider one spectral line, for 
which we assume, for simplicity, a Doppler profile centered at a 



frequency vq in the particle's frame. With these assumptions, the 
Doppler width Av D is a constant, and the monochromatic spe- 
cific intensity of the radiation emerging from the layer toward 
the observer is simply given by: 

I(v) = (l - r^W) Bv {T) (33) 

where rn(v) is the monochromatic optical depth from z\ to zi'- 

t«(v)= K{v,z)dz, (34) 
where the absorption coefficient k(v, z) is written as 

k(v,z) = Ke~^^) , (35) 

where v ctl (z), defined in Section 15.31 by Equation dT6b . is the 
shifted line-center frequency for a velocity V(z) at position z. 
Vctr(z) - v o (1 + V(z)/c), The quantity K is a constant that con- 
tains all the attributes of the line (population of the lower level, 
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Fig. 8. Comparison with ID plane-parallel models. The figures shows the monochromatic radiative quantities, at hv - 296 eV 
(A — 4.19 nm), calculated with IRIS (solid lines) and with a ID solver (symbols), in the case of a radiative shock structure made of 
homogeneous plane-parallel layers. The top left panel displays the specific intensity for different polar angles ft 0° (green, square), 
45° (blue, triangle), 89° (red, filled circle), 135° (orange, diamond), 180° (cyan, asterisk), and for a set of twelve azimuthal angles 
(p(°) = (0, 45, 89, 90, 135, 179, 180, 225, 269, 270, 315, 359). There is no legend for the <p values, because all the curves, for a given 
ft are perfectly superposed. The top right, bottom left, and bottom right panels show respectively the mean intensity / (black, 
square), the components of the radiative flux vector F (black, square), and the components of the radiation pressure tensor P 
(P zz : blue, square; P xx = P vy : red, plus). See Section [iT2"l for more details. 



and oscillator strength) and the value of the constant Doppler 
width. 

Now, since we use monotonic laws (cubic Hermite polyno- 
mials, see Section [572b to interpolate physical quantities along a 
short-characteristic, T\z(y) may be written as 



T«(V)= R'(V C 
•Jvar, 



-) K e \ Av ° I dv a 



(36) 



where R is the function relating the shifted transition frequency 
v ctr to the position z within the layer: z = R(v ca ), R' being its 
derivative, and where we define: 

|Vctr 2 = Vo[l + -f) 

In this simplified configuration, the optical depth appears as 
the convolution of a function A(v) and a Gaussian function B(v) 
(the line profile), both defined as follows: 

ti 2 (v) =(A*B) (y) , (38a) 

where 

R'(y ctr ) if Vctrj < v < v ctl - 2 



A(v) 







if v < y ctri or V > V, 



(38b) 



B(v) = Ke 



(38c) 



For the test, we adopt the following typical values of ve- 
locities in winds of early massive stars: V\ = 2000 kms~', 
V2 = 2250 km s _1 , and a thermal velocity (constant in our uni- 
form layer) of V t h = 25 km s~'. With such a choice, any line 
is Doppler shifted by an amount of 10 times its Doppler width 



between position zi and position zi, i.e., (v ctr , 



,)/Av D 



(V 2 - Vi)/Vth = 10. We also adopt a temperature T = 10 5 K, 
and a layer size Z2 — Zi = 10 9 cm. 

We consider an artificial line, arbitrarily centered at hvo = 
80.00 eV in the particle's frame. In the observer's frame, the line 



is centered at v c 



80.53 eV at zu and at v, 



Ctl"2 



80.60 eV at 



Ctl"2 



Z2- Note also that we have defined the K value (see Equationl35ll 
so that the transmission of the layer at the line center equals 0.5 
in the case when there is no velocity gradient. In this case in- 
deed, the monochromatic thermal absorption k(v, z) is indepen- 
dent of z, so that K is inferred from e~ K ^ Z2 ~ Zl ' - 0.5. This way, 
we avoid the extreme cases of an optically thin or an optically 
thick layer. The red thick curve in Figure [9] bottom left panel, 
displays the monochromatic specific intensity I(v) that emerges 
from the layer at position zi toward the observer, as calculated by 
IRIS. This emission spectrum is displayed versus photon energy 
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hv (eV) or photon wavelength A (A). For comparison, the black 
thick curve represents the line, centered at 80.53 eV, that would 
obtain if there was no velocity gradient, i.e., if V2 = Vi. The bot- 
tom right panel of the figure portrays the transmission, e~ Tl2(v) , 
of the layer. In the highly simplified configuration of our test, the 
transmission and the specific intensity are related by the linear 
Equation d3~3l >. This is why this transmission spectrum looks the 
same, but inverted, as the emission spectrum. In the presence of 
a velocity gradient, the absorption of the medium decreases, is 
shifted (blueward if the material approaches the observer, red- 
ward if it recedes from the observer), and spreads over a larger 
spectral range. This also applies to the emission of the medium. 

In this ideal test case, we can explain the shape and the po- 
sition of the line depicted by the red curve. The optical depth is 
given by Equation d36b . In addition, the layer is modeled with 
a unique cell. So, the cubic monotonic Hermite polynomial that 
defines V(z) along the short-characteristic between z\ and zi al- 
most coincides with a linear law, since the downwind endpoint 
is a ghost point in which state parameters are defined by linear 
extrapolations (see Section lBTTT ). In these specific conditions, z is 
trivially related to v ctr by a linear function R(v ca ), so that A(v) is 
merely a boxcar function. Then, the optical depth is 

c zi — 71 r Vc "2 J "ca-A 2 
Tia(v) = „ e\»») dv* (39) 



where v ctll and v ctr , are defined by Equations ( 1371 ). The calcula- 
tion of the monochromatic specific intensity with Equation ( T33l . 
using the relation ( 1391 above, provides a result that exactly fits 
the red curve obtained with IRIS. Note that this spectral line is 
centered at (v ctr] + v ctf2 )/2 » 80.57 eV. 

Last but not least, Fig. [9] (emission spectrum, left, or 
transmission spectrum, right) shows that the e D parameter (cf. 
Equation ( f20b ) plays a crucial role in the precision of the calcu- 
lated line, through the subgriding of the short-characteristic. The 
cyan curve shows the profile obtained with e D = 11, for which 
we checked that the short-characteristic is not subdivided. This 
profile is composed of two peaks, which are located in the two 
extreme positions of the line center, v ctri and v ct , 2 , with a large 
trough in the middle. Such a discrepancy between this profile 
and the correct red one may result in erroneous interpretations 
of observed spectra when compared with synthetic spectra. As 
e D decreases, the trough is being gradually filled, as shown by 
the other curves obtained with different e D : 6 (green), 5 (orange), 
3 (blue), 2 (grey). This test verifies that the line shape starts to 
stabilize from e D = 1. Tiny waves remain in the profile calcu- 
lated with e D — 1. This profile is not shown in the figure be- 
cause these waves are barely perceptible, so that it almost coin- 
cides with the red profile. Our numerical tests indicate that the 
short-characteristic is subdivided in 14 subintervals for e D — 1, 
whereas it is subdivided in 46 subintervals for e D = 0.3. Not 
surprisingly, increasing the precision involves a larger computa- 
tional cost. Such a test confirms that the ideal value of e D lies 
somewhere between 1/3 and 1 (see Section l53l l. 

We stress again that the profile depicted by the red curves in 
Fig.|9]has been obtained in the ideal case of a single line with a 
Doppler profile, under the assumption of a unique layer at uni- 
form temperature, in LTE, and with a linear velocity variation 
within the layer. In the general case, none of these situations 
holds, which results in a more complex profile. Only a full cal- 
culation with our subgriding method can provide the actual pro- 
file. 



9. Conclusion 

In this paper, we have described IRIS, a new generic three- 
dimensional spectral radiative transfer code, which solves the 
monochromatic static 3D radiative transfer equation in the ob- 
server's frame, in Cartesian coordinates. We drop the time 
derivative in the transfer equation because we consider situa- 
tions in which the dynamical timescales are large compared to 
the photon free-flight time, such that the studied medium is as- 
sumed to be in a frozen state during the time when the radia- 
tion field propagates throughout its structure. The code is pri- 
marily intended for post-processing any hydrodynamics snap- 
shot, i.e., any structure provided at a given instant by a (radi- 
ation)(magneto)hydrodynamics simulation. That way, IRIS can 
determine, at each instant, the radiation field of a non-stationary 
flow. We consider non-relativistic flow velocities. Currently, we 
assume local thermodynamic equilibrium. 

IRIS is mainly a diagnostic tool to be used for comparisons 
between predicted spectra, maps or images one the one hand, 
and astrophysical observations or laboratory astrophysics mea- 
surements on the other hand. It computes the monochromatic 
specific intensity and optical depth at any position within a hy- 
drodynamics structure, for any direction. Associated with appro- 
priate opacities specified by the user in a dedicated module, it 
calculates synthetic spectra and maps or images emerging from 
the studied structure or astrophysical object. It can also calcu- 
late, from the specific intensity and through angular integrations, 
the moments of the radiation field: the mean intensity, the radia- 
tion flux vector, and the radiation pressure tensor, at any position 
within the studied medium. 

IRIS works with any 3D Cartesian grid, the latter be- 
ing nonuniform in each direction. The code uses a short- 
characteristics solver to determine the formal solution of the ra- 
diative transfer equation. We have implemented a very efficient 
piecewise cubic, locally monotonic, interpolation technique, that 
considerably reduces the numerical diffusion effects of the short- 
characteristics method. The latter is used for interpolating any 
physical quantity in the cell faces of the computational grid, and 
for defining laws of variation of physical quantities along short- 
characteristics. 

IRIS is able to handle horizontal periodic boundary con- 
ditions of a simulation box. This configuration occurs for a 
medium with an infinite (that is, large compared to the extent 
of the computational box) extension in its horizontal plane and a 
double periodicity in this plane, along with a finite extension in 
its vertical direction, such as a stellar atmosphere or an accretion 
disk. 

Since it is formulated in the observer's frame, IRIS can deal 
with any (non-relativistic) non-monotonic macroscopic veloci- 
ties. The code can be applied to a large number of radiating as- 
trophysical objects or structures, such as accretion shocks or jets 
in young stellar objects, stellar atmospheres, exoplanet atmo- 
spheres, accretion disks, rotating stellar winds, and cosmolog- 
ical structures. IRIS has already been applied to predict X-UV 
spectra of la boratory generated radiative shocks (see conference 
proceedings Ibgui et al.ll20 1 2aUbt we are also writing a paper on 
this topic: Ibgui et al. 2012, in preparation). 

We envision various extensions in the near future. First, we 
will implement an iterative method based on the Accelerated 
Lambda Iteration (ALI) technique to handle scattering and 
NLTE effects. We will in terface IRIS with the NLTE stellar at - 
mosphere code TLUSTY (lHubenvll988HHubenv & Lanzll995l) . 
from which we take the treatment of atomic data, and all local 
physics (opacities, and the scheme for solving the set of statis- 
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Fig. 9. Test of the velocity gradient effect. The top panel represents a single plasma layer in a stellar wind, with a velocity gradient. 
The layer is at uniform temperature T, with one radiating species, assuming LTE and no scattering. The layer is modeled in IRIS 
with one cell between grid points zi and zz, with wind velocities V\ > (the material approaches the observer) and Vz > V\. The 
red thick curve in the bottom left panel displays the monochromatic specific intensity 7 V , due to the radiation that emerges from the 
layer at zz toward the observer, as computed by IRIS. This emission spectrum is plotted versus photon energy hv (eV) or photon 
wavelength A (A). It is obtained with a parameter e D = 0.3. Too large values of e D result in incorrect profiles, as shown by the 
following curves: cyan (e D = 11), green (e D = 6), orange (e D = 5), blue (e D = 3), grey (e D = 2). The bottom right panel displays 
the corresponding transmission spectrum. For comparison, the black thick curve portrays the spectral line that would obtain if there 
was no velocity gradient (V2 = Vi). See Section liOl for detailed explanations. 



tical equilibrium equations using the preconditioning method to 
obtain NLTE level populations). Later, extending IRIS to polar- 
ized radiative transfer will provide the code with the capability 
to diagnose magnetic fields. 

On the computational side, IRIS is written in Fortran 95. The 
current version runs on a single processor. We plan to paral- 
lelize the code with the message passing interface (MPI) library. 
Post-processing of (R)(M)HD calculations performed on adap- 
tively refined grids generated with the adaptive mesh refinement 
(AMR) technique will benefit from 3D nonuniform Cartesian 
grids in IRIS, though its full implementation will require addi- 
tional work. 
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